/*-
 * Copyright (c) 2013 Shang-Rong Cai. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 * $ MTBLAS TEST: tst_lv_3_gmm.c    Rev. 0.0.14    2013-11-11T12:25+08:00 $
 */

#include <stdlib.h>
#include <stdio.h>
#include <complex.h>

#include <cblas.h>

#include <tst_lv_a_com.h>

static const int flg_trn[3] = {MTBLAS_TRANS_NONE, MTBLAS_TRANS, MTBLAS_TRANS_CONJ};


int test_gemmf(size_t thrz)
{
    int ierr = ERRNO_UNDEFINED;
    int ii, ij;
    int cnt;

    float alpha, beta;
    float vc_alpha, vc_beta;

    vec_ff_t vx;
    vec_ff_t vy;
    vec_ff_t vz;

    mat_ff_t ma;
    mat_ff_t mb;
    mat_ff_t mc;
    mat_ff_t md;

    float buf_st;  /* Single-thread return */
    //float buf_mt;  /* Multi-threads return */
    float buf_cf;  /* cblas return */

    printf("gemmf() --> cblas_sgemm()\n");
    alpha = rand_coef_ff();
    beta = rand_coef_ff();
    printf("ALPHA = % .6E\n", alpha);
    printf("BETA = % .6E\n", beta);
    vc_alpha = 1.0;
    vc_beta = 1.0;

    ierr = mat_ff_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &ma);
    ierr = mat_ff_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mb);

    ierr = mat_ff_twins_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mc, &md);

    ierr = vec_ff_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_C_SIZE_RUN, &vx);
    ierr = vec_ff_twins_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_R_SIZE_RUN, &vy, &vz);

    cnt = 1;

    /* Row-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_ff_twins_reset(&mc, &md);
            ierr = vec_ff_twins_reset(&vy, &vz);

            ierr = stblas_gemmf(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = stblas_gemvf(MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_nrm2f(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_ff_twins_reset(&mc, &md);
            ierr = vec_ff_twins_reset(&vy, &vz);

            ierr = mtblas_gemmf(thrz, MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ii], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_gemvf(thrz, MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_nrm2f(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            ierr = mat_ff_twins_reset(&mc, &md);
            ierr = vec_ff_twins_reset(&vy, &vz);

            cblas_sgemm(CblasRowMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            cblas_sgemv(CblasRowMajor, CblasNoTrans, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            buf_cf = cblas_snrm2(PROPER_R_SIZE_RUN, vz.data, vz.step);

            printf("%2d. ROW_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_ff_stdout("    ST:", buf_st);
            //check_value_ff_stdout("    MT:", buf_mt);
            check_value_ff_stdout("    CF:", buf_cf);

            cnt++;
        }
    }

    /* Column-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_ff_twins_reset(&mc, &md);
            ierr = vec_ff_twins_reset(&vy, &vz);

            ierr = stblas_gemmf(MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            ierr = stblas_gemvf(MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_nrm2f(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_ff_twins_reset(&mc, &md);
            ierr = vec_ff_twins_reset(&vy, &vz);

            ierr = mtblas_gemmf(thrz, MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ii], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_gemvf(thrz, MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_nrm2f(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);

*/
            ierr = mat_ff_twins_reset(&mc, &md);
            ierr = vec_ff_twins_reset(&vy, &vz);

            cblas_sgemm(CblasColMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            cblas_sgemv(CblasColMajor, CblasNoTrans, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            buf_cf = cblas_snrm2(PROPER_R_SIZE_RUN, vz.data, vz.step);

            printf("%2d. COL_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_ff_stdout("    ST:", buf_st);
            //check_value_ff_stdout("    MT:", buf_mt);
            check_value_ff_stdout("    CF:", buf_cf);

            cnt++;
        }
    }

    ierr = vec_ff_delete(&vx);
    ierr = vec_ff_delete(&vy);
    ierr = vec_ff_delete(&vz);

    ierr = mat_ff_delete(&ma);
    ierr = mat_ff_delete(&mb);
    ierr = mat_ff_delete(&mc);
    ierr = mat_ff_delete(&md);

    return ierr;
}


int test_gemm(size_t thrz)
{
    int ierr = ERRNO_UNDEFINED;
    int ii, ij;
    int cnt;

    double alpha, beta;
    double vc_alpha, vc_beta;

    vec_df_t vx;
    vec_df_t vy;
    vec_df_t vz;

    mat_df_t ma;
    mat_df_t mb;
    mat_df_t mc;
    mat_df_t md;

    double buf_st;  /* Single-thread return */
    //double buf_mt;  /* Multi-threads return */
    double buf_cf;  /* cblas return */

    printf("gemm() --> cblas_dgemm()\n");
    alpha = rand_coef_df();
    beta = rand_coef_df();
    printf("ALPHA = % .15lE\n", alpha);
    printf("BETA = % .15lE\n", beta);
    vc_alpha = 1.0;
    vc_beta = 1.0;

    ierr = mat_df_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &ma);
    ierr = mat_df_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mb);

    ierr = mat_df_twins_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mc, &md);

    ierr = vec_df_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_C_SIZE_RUN, &vx);
    ierr = vec_df_twins_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_R_SIZE_RUN, &vy, &vz);

    cnt = 1;

    printf("Row-Major:\n");
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_df_twins_reset(&mc, &md);
            ierr = vec_df_twins_reset(&vy, &vz);

            ierr = stblas_gemm(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = stblas_gemv(MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_nrm2(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_df_twins_reset(&mc, &md);
            ierr = vec_df_twins_reset(&vy, &vz);

            ierr = mtblas_gemm(thrz, MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_gemv(thrz, MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_nrm2(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            ierr = mat_df_twins_reset(&mc, &md);
            ierr = vec_df_twins_reset(&vy, &vz);

            cblas_dgemm(CblasRowMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            cblas_dgemv(CblasRowMajor, CblasNoTrans, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            buf_cf = cblas_dnrm2(PROPER_R_SIZE_RUN, vz.data, vz.step);

            printf("%2d. ROW_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_df_stdout("    ST:", buf_st);
            //check_value_ff_stdout("    MT:", buf_mt);
            check_value_df_stdout("    CF:", buf_cf);

            cnt++;
        }
    }

    printf("Column-Major\n");
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_df_twins_reset(&mc, &md);
            ierr = vec_df_twins_reset(&vy, &vz);

            ierr = stblas_gemm(MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            ierr = stblas_gemv(MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_nrm2(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_df_twins_reset(&mc, &md);
            ierr = vec_df_twins_reset(&vy, &vz);

            ierr = mtblas_gemm(thrz, MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            ierr = mtblas_gemv(thrz, MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_nrm2(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            ierr = mat_df_twins_reset(&mc, &md);
            ierr = vec_df_twins_reset(&vy, &vz);

            cblas_dgemm(CblasColMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            cblas_dgemv(CblasColMajor, CblasNoTrans, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            buf_cf = cblas_dnrm2(PROPER_R_SIZE_RUN, vz.data, vz.step);

            printf("%2d. COL_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_df_stdout("    ST:", buf_st);
            //check_value_df_stdout("    MT:", buf_mt);
            check_value_df_stdout("    CF:", buf_cf);

            cnt++;
        }
    }

    ierr = vec_df_delete(&vx);
    ierr = vec_df_delete(&vy);
    ierr = vec_df_delete(&vz);

    ierr = mat_df_delete(&ma);
    ierr = mat_df_delete(&mb);
    ierr = mat_df_delete(&mc);
    ierr = mat_df_delete(&md);

    return ierr;
}


int test_gemml(size_t thrz)
{
    int ierr = ERRNO_UNDEFINED;
    int ii, ij;
    int cnt;

    long double alpha, beta;
    long double vc_alpha, vc_beta;

    vec_lf_t vx;
    vec_lf_t vy;
    vec_lf_t vz;

    mat_lf_t ma;
    mat_lf_t mb;
    mat_lf_t mc;
    mat_lf_t md;

    long double buf_st;  /* Single-thread return */
    //long double buf_mt;  /* Multi-threads return */

    printf("gemml()\n");
    alpha = rand_coef_lf();
    beta = rand_coef_lf();
    printf("ALPHA = % .24LE\n", alpha);
    printf("BETA = % .24LE\n", beta);
    vc_alpha = 1.0;
    vc_beta = 1.0;

    ierr = mat_lf_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &ma);
    ierr = mat_lf_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mb);

    ierr = mat_lf_twins_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mc, &md);

    ierr = vec_lf_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_C_SIZE_RUN, &vx);
    ierr = vec_lf_twins_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_R_SIZE_RUN, &vy, &vz);

    cnt = 1;

    /* Row-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_lf_twins_reset(&mc, &md);
            ierr = vec_lf_twins_reset(&vy, &vz);

            ierr = stblas_gemml(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = stblas_gemvl(MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_nrm2l(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_lf_twins_reset(&mc, &md);
            ierr = vec_lf_twins_reset(&vy, &vz);

            ierr = mtblas_gemml(thrz, MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_gemvl(thrz, MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_nrm2l(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            printf("%2d. ROW_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_lf_stdout("    ST:", buf_st);
            //check_value_lf_stdout("    MT:", buf_mt);

            cnt++;
        }
    }

    /* Column-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_lf_twins_reset(&mc, &md);
            ierr = vec_lf_twins_reset(&vy, &vz);

            ierr = stblas_gemml(MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            ierr = stblas_gemvl(MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_nrm2l(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_lf_twins_reset(&mc, &md);
            ierr = vec_lf_twins_reset(&vy, &vz);

            ierr = mtblas_gemmf(thrz, MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            ierr = mtblas_gemvf(thrz, MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_nrm2f(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            printf("%2d. COL_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_lf_stdout("    ST:", buf_st);
            //check_value_lf_stdout("    MT:", buf_mt);

            cnt++;
        }
    }

    ierr = vec_lf_delete(&vx);
    ierr = vec_lf_delete(&vy);
    ierr = vec_lf_delete(&vz);

    ierr = mat_lf_delete(&ma);
    ierr = mat_lf_delete(&mb);
    ierr = mat_lf_delete(&mc);
    ierr = mat_lf_delete(&md);

    return ierr;
}


int test_cgemmf(size_t thrz)
{
    int ierr = ERRNO_UNDEFINED;
    int ii, ij;
    int cnt;

    float complex alpha, beta;
    float complex vc_alpha, vc_beta;

    vec_fc_t vx;
    vec_fc_t vy;
    vec_fc_t vz;

    mat_fc_t ma;
    mat_fc_t mb;
    mat_fc_t mc;
    mat_fc_t md;

    float buf_st;  /* Single-thread return */
    //float complex buf_mt;  /* Multi-threads return */
    float buf_cf;  /* cblas return */

    printf("cgemmf() --> cblas_cgemm()\n");
    alpha = rand_coef_fc();
    beta = rand_coef_fc();
    printf("ALPHA = % .6E + i % .6E\n", crealf(alpha), cimagf(alpha));
    printf("BETA = % .6E + i % .6E\n", crealf(beta), cimagf(beta));
    vc_alpha = 1.0;
    vc_beta = 1.0;

    ierr = mat_fc_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &ma);
    ierr = mat_fc_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mb);

    ierr = mat_fc_twins_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mc, &md);

    ierr = vec_fc_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_C_SIZE_RUN, &vx);
    ierr = vec_fc_twins_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_R_SIZE_RUN, &vy, &vz);

    cnt = 1;

    /* Row-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_fc_twins_reset(&mc, &md);
            ierr = vec_fc_twins_reset(&vy, &vz);

            ierr = stblas_cgemmf(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = stblas_cgemvf(MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_cnrm2f(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_fc_twins_reset(&mc, &md);
            ierr = vec_fc_twins_reset(&vy, &vz);

            ierr = mtblas_cgemmf(thrz, MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_cgemvf(thrz, MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_cnrm2f(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            ierr = mat_fc_twins_reset(&mc, &md);
            ierr = vec_fc_twins_reset(&vy, &vz);

            cblas_cgemm(CblasRowMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, (void*) &alpha, ma.data, ma.row_size, mb.data, mb.row_size, (void*) &beta, md.data, md.row_size);
            cblas_cgemv(CblasRowMajor, CblasNoTrans, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, (void*) &vc_alpha, md.data, md.row_size, vx.data, vx.step, (void*) &vc_beta, vz.data, vz.step);
            buf_cf = cblas_scnrm2(PROPER_R_SIZE_RUN, vz.data, vz.step);

            printf("%2d. ROW_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_ff_stdout("    ST:", buf_st);
            //check_value_ff_stdout("    MT:", buf_mt);
            check_value_ff_stdout("    CF:", buf_cf);

            cnt++;
        }
    }

    /* Column-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_fc_twins_reset(&mc, &md);
            ierr = vec_fc_twins_reset(&vy, &vz);

            ierr = stblas_cgemmf(MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            ierr = stblas_cgemvf(MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_cnrm2f(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_fc_twins_reset(&mc, &md);
            ierr = vec_fc_twins_reset(&vy, &vz);

            ierr = mtblas_cgemmf(thrz, MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_cgemvf(thrz, MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_cnrm2f(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            ierr = mat_fc_twins_reset(&mc, &md);
            ierr = vec_fc_twins_reset(&vy, &vz);

            cblas_cgemm(CblasColMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, (void*) &alpha, ma.data, ma.col_size, mb.data, mb.col_size, (void*) &beta, md.data, md.col_size);
            cblas_cgemv(CblasColMajor, CblasNoTrans, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, (void*) &vc_alpha, md.data, md.col_size, vx.data, vx.step, (void*) &vc_beta, vz.data, vz.step);
            buf_cf = cblas_scnrm2(PROPER_R_SIZE_RUN, vz.data, vz.step);

            printf("%2d. COL_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_ff_stdout("    ST:", buf_st);
            //check_value_ff_stdout("    MT:", buf_mt);
            check_value_ff_stdout("    CF:", buf_cf);

            cnt++;
        }
    }

    ierr = vec_fc_delete(&vx);
    ierr = vec_fc_delete(&vy);
    ierr = vec_fc_delete(&vz);

    ierr = mat_fc_delete(&ma);
    ierr = mat_fc_delete(&mb);
    ierr = mat_fc_delete(&mc);
    ierr = mat_fc_delete(&md);

    return ierr;
}


int test_cgemm(size_t thrz)
{
    int ierr = ERRNO_UNDEFINED;
    int ii, ij;
    int cnt;

    double complex alpha, beta;
    double complex vc_alpha, vc_beta;

    vec_dc_t vx;
    vec_dc_t vy;
    vec_dc_t vz;

    mat_dc_t ma;
    mat_dc_t mb;
    mat_dc_t mc;
    mat_dc_t md;

    double buf_st;  /* Single-thread return */
    //double complex buf_mt;  /* Multi-threads return */
    double buf_cf;  /* cblas return */

    printf("cgemm() --> cblas_zgemm()\n");
    alpha = random_double() - random_double() + _Complex_I * (random_double() - random_double());
    beta = random_double() - random_double() + _Complex_I * (random_double() - random_double());
    printf("ALPHA = % .15lE + i % .15lE\n", creal(alpha), cimag(alpha));
    printf("BETA = % .15lE + i % .15lE\n", creal(beta), cimag(beta));
    vc_alpha = 1.0;
    vc_beta = 1.0;

    ierr = mat_dc_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &ma);
    ierr = mat_dc_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mb);

    ierr = mat_dc_twins_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mc, &md);

    ierr = vec_dc_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_C_SIZE_RUN, &vx);
    ierr = vec_dc_twins_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_R_SIZE_RUN, &vy, &vz);

    cnt = 1;

    /* Row-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_dc_twins_reset(&mc, &md);
            ierr = vec_dc_twins_reset(&vy, &vz);

            ierr = stblas_cgemm(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = stblas_cgemv(MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_cnrm2(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_dc_twins_reset(&mc, &md);
            ierr = vec_dc_twins_reset(&vy, &vz);

            ierr = mtblas_cgemm(thrz, MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_cgemv(thrz, MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_cnrm2(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            ierr = mat_dc_twins_reset(&mc, &md);
            ierr = vec_dc_twins_reset(&vy, &vz);

            cblas_zgemm(CblasRowMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, (void*) &alpha, ma.data, ma.row_size, mb.data, mb.row_size, (void*) &beta, md.data, md.row_size);
            cblas_zgemv(CblasRowMajor, CblasNoTrans, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, (void*) &vc_alpha, md.data, md.row_size, vx.data, vx.step, (void*) &vc_beta, vz.data, vz.step);
            buf_cf = cblas_dznrm2(PROPER_R_SIZE_RUN, vz.data, vz.step);

            printf("%2d. ROW_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_df_stdout("    ST:", buf_st);
            //check_value_df_stdout("    MT:", buf_mt);
            check_value_df_stdout("    CF:", buf_cf);

            cnt++;
        }
    }

    /* Column-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_dc_twins_reset(&mc, &md);
            ierr = vec_dc_twins_reset(&vy, &vz);

            ierr = stblas_cgemm(MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            ierr = stblas_cgemv(MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_cnrm2(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_dc_twins_reset(&mc, &md);
            ierr = vec_dc_twins_reset(&vy, &vz);

            ierr = mtblas_cgemm(thrz, MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_cgemv(thrz, MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_cnrm2(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            ierr = mat_dc_twins_reset(&mc, &md);
            ierr = vec_dc_twins_reset(&vy, &vz);

            cblas_zgemm(CblasColMajor, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, (void*) &alpha, ma.data, ma.col_size, mb.data, mb.col_size, (void*) &beta, md.data, md.col_size);
            cblas_zgemv(CblasColMajor, CblasNoTrans, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, (void*) &vc_alpha, md.data, md.col_size, vx.data, vx.step, (void*) &vc_beta, vz.data, vz.step);
            buf_cf = cblas_dznrm2(PROPER_R_SIZE_RUN, vz.data, vz.step);

            printf("%2d. COL_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_df_stdout("    ST:", buf_st);
            //check_value_df_stdout("    MT:", buf_mt);
            check_value_df_stdout("    CF:", buf_cf);

            cnt++;
        }
    }

    ierr = vec_dc_delete(&vx);
    ierr = vec_dc_delete(&vy);
    ierr = vec_dc_delete(&vz);

    ierr = mat_dc_delete(&ma);
    ierr = mat_dc_delete(&mb);
    ierr = mat_dc_delete(&mc);
    ierr = mat_dc_delete(&md);

    return ierr;
}

int test_cgemml(size_t thrz)
{
    int ierr = ERRNO_UNDEFINED;
    int ii, ij;
    int cnt;

    long double complex alpha, beta;
    long double complex vc_alpha, vc_beta;

    vec_lc_t vx;
    vec_lc_t vy;
    vec_lc_t vz;

    mat_lc_t ma;
    mat_lc_t mb;
    mat_lc_t mc;
    mat_lc_t md;

    long double buf_st;  /* Single-thread return */
    //double complex buf_mt;  /* Multi-threads return */

    printf("cgemml()\n");
    alpha = random_long_double() - random_long_double() + _Complex_I * (random_long_double() - random_long_double());
    beta = random_long_double() - random_long_double() + _Complex_I * (random_long_double() - random_long_double());
    printf("ALPHA = % .24LE + i % .24LE\n", creall(alpha), cimagl(alpha));
    printf("BETA = % .24LE + i % .24LE\n", creall(beta), cimagl(beta));

    vc_alpha = 1.0;
    vc_beta = 1.0;

    ierr = mat_lc_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &ma);
    ierr = mat_lc_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mb);

    ierr = mat_lc_twins_create(PROPER_MATRIX_MAX, PROPER_MATRIX_MIN, &mc, &md);

    ierr = vec_lc_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_C_SIZE_RUN, &vx);
    ierr = vec_lc_twins_create(PROPER_VECTOR_MAX, PROPER_VECTOR_MIN, PROPER_R_SIZE_RUN, &vy, &vz);

    cnt = 1;

    /* Row-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_lc_twins_reset(&mc, &md);
            ierr = vec_lc_twins_reset(&vy, &vz);

            ierr = stblas_cgemml(MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = stblas_cgemvl(MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_cnrm2l(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_lc_twins_reset(&mc, &md);
            ierr = vec_lc_twins_reset(&vy, &vz);

            ierr = mtblas_cgemml(thrz, MTBLAS_ROW_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_cgemvl(thrz, MTBLAS_ROW_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_cnrm2l(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            printf("%2d. ROW_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_lf_stdout("    ST:", buf_st);
            //check_value_lf_stdout("    MT:", buf_mt);

            cnt++;
        }
    }

    /* Column-Major */
    for (ii = 0; ii < 3; ii++)
    {
        for (ij = 0; ij < 3; ij++)
        {
            ierr = mat_lc_twins_reset(&mc, &md);
            ierr = vec_lc_twins_reset(&vy, &vz);

            ierr = stblas_cgemml(MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.col_size, mb.data, mb.col_size, beta, md.data, md.col_size);
            ierr = stblas_cgemvl(MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.col_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = stblas_cnrm2l(PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
/*
            ierr = mat_lc_twins_reset(&mc, &md);
            ierr = vec_lc_twins_reset(&vy, &vz);

            ierr = mtblas_cgemml(thrz, MTBLAS_COL_MAJOR, flg_trn[ii], flg_trn[ij], PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, PROPER_K_SIZE_RUN, alpha, ma.data, ma.row_size, mb.data, mb.row_size, beta, md.data, md.row_size);
            ierr = mtblas_cgemvl(thrz, MTBLAS_COL_MAJOR, MTBLAS_TRANS_NONE, PROPER_R_SIZE_RUN, PROPER_C_SIZE_RUN, vc_alpha, md.data, md.row_size, vx.data, vx.step, vc_beta, vz.data, vz.step);
            ierr = mtblas_cnrm2l(thrz, PROPER_R_SIZE_RUN, vz.data, vz.step, &buf_st);
*/
            printf("%2d. COL_MAJOR, TRANSPOSE FLAG: %d, TRANSPOSE FLAG: %d\n", cnt, ii, ij);
            check_value_lf_stdout("    ST:", buf_st);
            //check_value_lf_stdout("    MT:", buf_mt);

            cnt++;
        }
    }

    ierr = vec_lc_delete(&vx);
    ierr = vec_lc_delete(&vy);
    ierr = vec_lc_delete(&vz);

    ierr = mat_lc_delete(&ma);
    ierr = mat_lc_delete(&mb);
    ierr = mat_lc_delete(&mc);
    ierr = mat_lc_delete(&md);

    return ierr;
}
